A  Two-Dimensional  Meteorological  Computer  Model  for  the 

Forest  Canopy 


by  Arnold  Tunick 


ARL-MR-569 


August  2003 


Approved  for  public  release;  distribution  unlimited. 


NOTICES 

Disclaimers 

The  findings  in  this  report  are  not  to  be  construed  as  an  official 
Department  of  the  Army  position,  unless  so  designated  by  other 
authorized  documents. 

Citation  of  manufacturers’  or  trade  names  does  not  constitute  an  official 
endorsement  or  approval  of  the  use  thereof. 


Army  Research  Laboratory 

Adelphi,  MD  20783-1197 


ARL-MR-569 _ August  2003 


A  Two-Dimensional  Meteorological  Computer  Model  for  the 

Forest  Canopy 

Arnold  Tunick 

Computational  and  Information  Sciences  Directorate,  ARL 


Approved  for  public  release;  distribution  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to 
comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

August  2003  Final 


4.  TITLE  AND  SUBTITLE 

A  Two-Dimensional  Meteorological  Computer  Model  for  the  Forest  Canopy 


3.  DATES  COVERED  (From  -  To) 

October  2002- June  2003 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Research  Laboratory 
Attn:  AMSRL-CI-EE 
2800  Powder  Mill  Road 
Adelphi,  MD  20783-1197 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Research  Laboratory 
2800  Powder  Mill  Road 
Adelphi,  MD  20783-1197 


5c.  PROGRAM  ELEMENT  NUMBER 

61102A 


5d,  PROJECT  NUMBER 

3FEJ00 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


ARL-MR-569 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 


13.  SUPPLEMENTARY  NOTES 

AMS  Code:  61 1 10253A1 1 
DA  Project:  B53A 


14.  ABSTRACT 

This  report  presents  the  equation  set,  modeling  assumptions,  and  some  initial  results  from  a  new,  physics-based  computer  model  that  is  being 
developed  for  two-dimensional  forest  canopy  wind  flow,  temperature,  and  turbulence  calculations.  The  model  is  based  on  the  conservation 
(simplified  Navier-Stokes)  equations  for  continuity,  momentum,  Reynolds  stress,  energy,  heat  flux,  and  turbulent  temperature  variance.  A  set 
of  simultaneous  equations  for  each  of  12  computed  variables  is  solved  iteratively  on  a  computational  grid  consisting  of  10  x  60  points. 
Horizontal  grid  spacing  is  50  m,  and  vertical  grid  spacing  is  0.5  m.  The  model  domain  is  500  x  30  m.  It  is  anticipated  that  improved 
turbulence  and  micrometeorological  models  for  forest  canopies  will  become  increasingly  useful  for  military  acoustic  application  research. 


15.  SUBJECT  TERMS 

Conservation  equations,  wind  flow,  energy  budget,  and  computational  methods 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Arnold  Tunick 

a.  REPORT 

UNCLASSIFIED 

b.  ABSTRACT 

UNCLASSIFIED 

c.  THIS  PAGE 

UNCLASSIFIED 

UL 

30 

19b.  TELEPHONE  NUMBER  ( Include  area  code) 
(301)394-1765 

Standard  Form  298  (Rev.  8/98) 

Prescribed  by  ANSI  Std.  Z39.18 


Contents 


Contents  iii 

List  of  Figures  iv 

Acknowledgments  v 

1.  Introduction  1 

2.  Forest  Canopy  Model  1 

2 . 1  Conservation  Equations . 1 

2.2  Modeling  Assumptions . 3 

2.3  Second-Order  Turbulence  Closure  Model  for  2-D  Forest  Canopies . 5 

2.4  Numerical  Methods . 8 

2.5  Forest  Canopy  Architecture . 10 

3.  Model  Results  10 

3 . 1  Uniform  Forest  Stands . 10 

3.2  Nonunifonn  Forest  Stands . 13 

4.  Summary  and  Conclusions  18 

5.  References  19 


iii 


List  of  Figures 


Figure  1.  Normalized  vertical  profiles  of  leaf  area  density  for  forest  canopies . 10 

Figure  2.  Model  results  for  uniform  forest  stands:  (a)  horizontal  wind  velocity,  (tij ,  in  units 
ms-1,  and  (b)  air  temperature,  ll)) ,  in  units  °C.  For  this  example,  canopy  height  ( h )  is 
10  m . 11 

Figure  3.  Profiles  of  horizontal  wind  velocity  ( (uj ,  in  units  ms-1)  and  air  temperature  ( (d\ ,  in 
units  °C)  derived  from  the  current  2-D  calculation  shown  in  comparison  to  1-D  profiles 
derived  from  an  earlier  study . 12 


Figure  4.  Spurious  oscillations  due  to  numerical  instability  in  the  2-D  model  results  for  uniform 
forest  stands.  The  computed  variable  is  horizontal  wind  velocity,  (u\ ,  in  units  ms~ 1 . 12 

Figure  5.  Model  results  for  nonuniform  forest  stands,  i.e.,  those  that  contain  a  single  step  change 
in  canopy  height:  (a)  horizontal  wind  velocity,  M  ,  in  units  ms'1,  and  (b)  wind  flow 
streamlines.  A  single  step  change  in  canopy  height  at  X/2  +  Ax  is  shown  using  open 
rectangles,  where  h  =  8  m  on  the  left  side  and  h  =  10  m  on  the  right  side . 13 

Figure  6.  Model  results  for  nonunifonn  forest  stands:  (a)  vertical  wind  velocity,  (wj ,  in  units 
ms'1  and  (b)  kinematic  (fluctuation)  pressure  ,  in  units  m  2s”2.  A  single  step  change  in 
canopy  height  occurs  at  X/2  +  Ax,  where  h  =  8  m  on  the  left  side  and  h  =  10  m  on  the  right 


side . 14 

Figure  7.  Model  results  for  nonunifonn  forest  stands:  (a)  air  temperature,  (o'j ,  in  units  °C  and 
(b)  effective  speed  of  sound,  (ceff  'j ,  in  units  ms'1.  A  single  step  change  in  canopy  height 
occurs  at  X/2  +  Ax,  where  h  =  8  m  on  the  left  side  and  h=  10  m  on  the  right  side . 14 

Figure  8.  Spurious  oscillations  due  to  numerical  instability  in  the  2-D  model  results  for 

nonunifonn  forest  stands.  The  computed  variable  is  horizontal  wind  velocity,  (u\  ,  in  units 
ms  1 . 15 

Figure  9.  Model  results  for  nonunifonn  forest  stands,  i.e.,  those  that  contain  a  single  step  change 
in  canopy  height:  (a)  horizontal  wind  velocity,  M  ,  in  units  ms'1,  and  (b)  wind  flow 
streamlines.  A  single  step  change  in  canopy  height  at  X/2  +  Ax  is  shown  using  open 
rectangles,  where  /i=4m  on  the  left  side  and  h  =  10  m  on  the  right  side . 16 

Figure  10.  Model  results  for  nonunifonn  forest  stands:  (a)  vertical  wind  velocity,  (zv\ ,  in  units 
ms'1,  and  (b)  kinematic  (fluctuation)  pressure  (j^ ,  in  units  m  2s'2.  A  single  step  change  in 
canopy  height  occurs  at  X/2  +  Ax,  where  h  =  4  m  on  the  left  side  and  h  =  10  m  on  the  right 
side . 16 


Figure  1 1 .  Model  results  for  nonunifonn  forest  stands:  (a)  air  temperature,  (d'j ,  in  units  °C,  and 
(b)  effective  speed  of  sound,  (ceff  'j ,  in  units  ms  1 .  A  single  step  change  in  canopy  height 
occurs  at  X/2  +  Ax,  where  h  =  4  m  on  the  left  side  and  h  =  10  m  on  the  right  side . 17 

Figure  12.  Numerical  instability,  i.e.,  2Ax  waves,  in  the  computation  of  wind  flow  streamlines 
for  nonunifonn  forest  stands:  (a)  jdiy  =  -[wdx  +  ^udz  and  (b)  j d <//  =  - j wdx  only.  A  single 
step  change  in  canopy  height  occurs  at  X/2  +  Ax,  where  h  =  4  m  on  the  left  side  and  h  =  10  m 
on  the  right  side . 17 


IV 


Acknowledgments 


The  author  would  like  to  thank  Ronald  Meyers,  Keith  Deacon,  and  Yansen  Wang  of  the  U.S. 
Army  Research  Laboratory  for  many  helpful  discussions  on  computational  methods.  The 
Department  of  the  Army  funded  this  research  through  Project  B53A/61 102A. 


v 


Intentionally  Left  Blank. 


vi 


1.  Introduction 


In  a  previous  study  ( 1 ,  2),  it  was  found  that  a  useful  mathematical  representation  of  the  wind 
flow,  temperatures,  and  turbulence  inside  and  above  a  (uniform)  continuous  forest  stand  could  be 
obtained  by  means  of  a  one-dimensional  (1-D),  steady-state,  second-order  turbulence  closure 
model,  with  an  embedded  radiative  transfer  and  energy  budget  algorithm  to  predict  the  heat 
source.  Development  of  this  model  made  it  possible  to  generate  realistic  profiles  for  effective 
sound  speed  inside  and  above  a  forest  canopy.  In  turn,  these  data  were  used  as  input  to  an 
acoustic  propagation  model  that  predicts  atmospheric  and  terrain  effects  on  short-range  acoustic 
attenuation  ( 3 ,  4).  Asa  result,  it  was  shown  that  attenuation  and  “ducting”  of  acoustic  waves  in 
and  around  forests  is  significantly  influenced  by  local  micrometeorological  profile  structure. 

However,  forest  stands  are  typically  inhomogeneous,  containing  nonunifonn  distributions  of 
canopy  height  and  leaf  area  density  (5).  In  addition,  open  fields,  roadways,  and  buildings  often 
border  forests.  Hence,  to  begin  to  address  nonunifonn  forests  and  forest  edges,  this  report 
presents  the  equation  set,  modeling  assumptions,  and  some  initial  results  from  a  new,  physics- 
based  computer  model  that  is  being  developed  for  two-dimensional  (2-D)  forest  canopy  wind 
flow,  temperature,  and  turbulence  calculations.  Like  the  earlier  1-D  model,  the  2-D  model  is 
based  on  the  conservation  (simplified  Navier-Stokes)  equations  for  continuity,  momentum, 
Reynolds  stress,  energy,  heat  flux,  and  turbulent  temperature  variance.  However,  in  this  case,  a 
set  of  simultaneous  equations  for  each  of  12  computed  variables  is  solved  iteratively  on  a 
computational  grid  consisting  of  10  x  60  points.  Horizontal  grid  spacing  is  50  m  and  vertical 
grid  spacing  is  0.5  m.  The  model  domain  is  500  x  30  m.  It  is  anticipated  that  improved  physics- 
based  theory  and  computer  modeling  for  meteorology  coupled  to  acoustics  will  become 
increasingly  useful  to  predict  effective  sound  speed  information  for  military  acoustic  application 
research  ( 6 ,  7). 

The  mathematical  model  for  the  forest  canopy  is  described  in  section  2.  Initial  model  results  for 
uniform  forest  stands  are  shown  in  section  3.1.  Initial  model  results  for  nonunifonn  forest 
stands,  i.e.,  those  that  contain  a  single  step  change  in  canopy  height,  are  shown  in  section  3.2.  A 
summary  and  conclusions  are  provided  in  section  4. 


2.  Forest  Canopy  Model 


2.1  Conservation  Equations 

The  conservation  (simplified  Navier  Stokes)  equations  for  the  current  model,  neglecting  coriolis* 
forces,  can  be  expressed  as  follows: 


Other  than  a  few  authors  (8-10),  most  consider  the  effect  of  the  coriolis  force  as  being  negligible  for  the  scales  of  motion 
considered. 
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The  continuity  equation  is 
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Here,  t  is  the  independent  variable  time,  p  is  air  density,  ui  is  the  i-component  of  the  wind 
velocity  vector  ( w,  =  u,  u2  =  v,  u-,  =  w),  and  x,  is  the  i-component  of  the  position  vector  (x\  =  x, 

xi  =y,  X3  =  z).  In  addition,/;  is  static  pressure  normalized  by  air  density  (i.e.,  kinematic 
pressure),  g  is  the  acceleration  due  to  gravity,  T  is  the  absolute  temperature,  do  is  a  deviation 
from  a  reference  temperature  that  decreases  with  height  at  the  adiabatic  lapse  rate,  and  v  is 
kinematic  viscosity.  The  overbar  and  primed  variables  indicate  the  mean  (time  averaged)  and 
fluctuating  components  of  the  given  quantity,  whereas  the  brackets,  ^  'j ,  and  double  primed 

variables  indicate  horizontal  averaging  and  departures  from  the  horizontal  averaging  operator 

(11). 


The  stress  equation  for  non-adiabatic  conditions  is 
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The  energy  equation  is 
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Here,  <9  is  ambient  air  temperature,  Kj  is  thermal  diffusivity,  and  So  is  the  forest  canopy  heat 
source  (or  sink).  The  heat  source  can  be  expressed  as  a  function  of  leaf  surface-to-ambient-air 
temperature  differences  (12). 
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The  heat  flux  equation  is 
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Finally,  the  turbulent  temperature  variance  equation  is 
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(6) 


2.2  Modeling  Assumptions 

Wilson  and  Shaw  (13)  give  the  following  closure  approximation  for  the  pressure  drag  force  in 
equation  2: 


%T~  =  C.A'flu,)  .  (7) 

o  xt  11 

2  —3 

Here,  Cd  is  the  forest  canopy  drag  coefficient  (=  0.10 )  and  A  (in  units  mm  )  is  the  leaf  area 
density.  It  is  further  assumed  that  pressure  forces  are  the  main  contributor  to  the  total  drag  from 
the  forest  canopy  (i.e.,  viscous  drag  forces  are  neglected). 

In  equation  3,  Mellor  (14)  and  Mellor  and  Yamada  (15,  16)  parameterize  the  triple-velocity 
products  as 
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where  the  turbulent  kinetic  energy  (t.k.e.)  is  q-  (iftift  ''j  and  /fi  is  a  function  of  mixing  length.* 

The  pressure-velocity  gradient  terms  (i.e.,  the  pressure  redistribution  terms)  in  equation  3  can  be 
rewritten  as 


(9) 


where 


*In  these  equations,  A i  through  /t7  are  length  scales,  which  contain  a  set  of  seven  closure  constants,  i.e.,  Ak  =  akl,  where  k  is  an 
arbitrary  index  1-7  and  /  is  the  mixing  length.  Values  for  these  closure  constants  are  as  follows:  at  =  0.39,  a2  =  0.85,  a}  =  16.57, 
=  0.23,  a5  =  0.74,  and  a2  =  10.10  (11.  13-16). 
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This  is  modeled  according  to  the  retum-to-isotropy  principle,  as  described  by  Mellor  {14)  and 
Donaldson  (17).  In  equation  10,  C  is  a  constant  whose  value  is  about  0.077  {11). 


Viscous  dissipation  is  assumed  isotropic  and  a  function  of  local  t.k.e.  intensity.  Viscous 
dissipation,  as  described  by  Katul  and  Albertson  {11),  is  parameterized  as 
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In  equation  4,  the  heat  source  term  (S'#),  as  described  by  Meyers  and  Paw  U  {12),  is  modeled  as 

S0=2A  i^-0)/rh  .  (13) 

Here,  A  (in  units  m2m-3)  is  the  leaf  area  density,  (o,  -  0  )  is  the  mean  leaf  surface-to-ambient-air 

temperature  difference,  and  rk  is  the  aerodynamic  resistance  to  heat  transfer.  A  1-D  radiative 
transfer  and  energy  budget  algorithm  is  incorporated  into  the  2-D  model  calculation  to  make  it 
possible  to  determine  the  heat  source  for  any  time  of  day.  To  do  this,  the  formulations  outlined 
by  Rachele  and  Tunick  {18)  are  used  to  calculate  the  incoming  total  radiation  at  the  canopy  top 
as  a  function  of  latitude,  longitude,  day  of  year,  and  time  of  day  (i.e.,  these  input  are  needed  to 
determine  the  solar  declination,  hour,  and  zenith  angles).  Then,  the  equations  provided  by  Weiss 
and  Nonnan  {19)  are  used  to  calculate  the  spectral  components  for  short-wave  (i.e.,  direct  beam 
and  diffuse,  visible  and  near-infrared)  radiation  as  a  function  of  the  total  downward  short-wave 
flux  at  canopy  top  because  extinction  and  reflection  through  the  forest  canopy  are  different  for 
each.  The  remainder  of  the  1-D  radiative  transfer  subroutine  for  the  forest  canopy  (i.e., 
transmission,  reflection,  absorption,  and  emission  of  the  solar  flux)  is  derived  from  the 
formulations  given  in  the  texts  by  Campbell  {20)  and  Campbell  and  Norman  {21). 

Returning  now  to  equations  5  and  6,  the  triple-velocity-temperature  products  contained  therein 
can  be  expressed  in  the  fonn  described  by  Mellor  {14)  and  Mellor  and  Yamada  {15,  16)  as 
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In  addition,  the  pressure  interaction  term  can  be  modeled  as 
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Finally,  the  molecular  dissipation  of  heat  can  be  modeled  as 
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This  completes  the  basic  equation  set  and  modeling  assumptions  for  the  2-D  forest  canopy 
model. 


2.3  Second-Order  Turbulence  Closure  Model  for  2-D  Forest  Canopies 

The  2-D  computer  model  currently  being  developed  is  relatively  unique.  This  is  because  higher- 
order  closure  models  reported  in  the  literature  (e.g.,  11-13,  22)  have  generally  focused  on 
estimates  of  1-D,  adiabatic  wind  flow,  and  turbulence  within  and  above  homogeneous  (unifonn) 
forest  canopies.  In  contrast,  the  new  model  may  be  applied  night  or  day  to  both  uniform  and 
nonuniform  stands  (and  forest  edges,  possibly).  The  new  model,  which  is  based  on  the  earlier 
works  of  Katul  and  Albertson  (11),  Meyers  and  Paw  U  (12),  Wilson  and  Shaw  (13),  and 
Donaldson  (17),  calculates  the  2-D,  steady  state,  canopy  wind  flow,  temperatures,  turbulent 
variances,  Reynolds  stress,  and  heat  flux.  The  parameterized  2-D  model  equations  for 
continuity,  the  mean  flow-longitudinal  (u) ,  the  mean  flow-vertical  (w) ,  Reynolds  stress 

( u'w ') ,  longitudinal  (i(2^ ,  lateral  (v'2^j ,  and  vertical  velocity  (u>'2^  variances,  the  mean 

temperature  ,  vertical  heat  flux  (iv'0'^j  ,  horizontal  heat  flux  (u' 0'^j ,  and  the  turbulent 

temperature  variance  /<9,2\ ,  respectively,  are  as  follows: 
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dx 

did'2 


qA6- 


dx 


—  { w 


dz 


d(0'  1  8( • - \  8(6 )  /— — r\  d(6 


-2 tu  6  H-^-2 (w  6 
dz  \  /  dx  \  I  dz 


qAfr- 


d(0'- 


8z 


Iff2 


(28) 


In  the  model,  mean  pressure  is  assumed  approximately  hydrostatic.  In  contrast,  the  kinematic 
(fluctuation)  pressure  l^'j  is  determined,  as  discussed  in  the  text  by  Ferziger  and  Peric  (23),  by 
taking  the  divergence  of  the  mean  flow  equations,  i.e., 


d 


dx 


V 


dt 


d_ 

dx 

V  J 


d(u'w' 

dz 


d\p_ 

dx2 


-CdA\u\- 


d(u 


dx 


(29) 


and 


d_ 

dz 


rd(wX ^ 


dt 


8z\  dx 


d(w ) 


w- 
dz  \  dz 


r 

d{u  w 


dx 


^)+^-CdA\u\- 


dz 


+  tLm 

T  dz 


dlw) 


dz 


(30) 


which,  after  some  rearranging  and  cancellation  of  tenns  due  to  continuity,  yields  Poisson’s 
equation  (i.e.,  pressure-velocity  coupling)  as 


d_ 

dx 

V  J 


d(u'w'\ 

dz 


d_\  d(u'w') 

dz  dx 
v  y 


d2("A)  ,  £  A) 

Sz 2  T  dz 


(31) 


Also,  note  that  the  lateral  velocity  variance  /v'2 )  is  calculated  in  the  current  2-D  model 

(equation  23)  to  support  the  turbulence  (t.k.e.)  closure  approximations  for  the  triple  product  and 
dissipation  terms. 
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2.4  Numerical  Methods 


A  computational  grid  consisting  of  10  x  60  points  is  chosen,  where  the  horizontal  grid  spacing 
( Ax  )  is  50  m  and  the  vertical  grid  spacing  ( Az )  is  0.5  m.  Hence,  the  model  domain  is  500  x 
30  m.  Vertical  derivatives  are  solved  using  a  lower-order  central  differencing  scheme  (23).  The 
first  derivatives  can  be  expressed  as 


8<Kuj) 

8z  2Az 


(32) 


and  the  second  derivatives  as 


8%,j)  ^j-i)-2^(/j)  +  ^(»j+i) 

0z2  (Az)2 


(33) 


Here,  i  is  the  horizontal  grid  index  and  j  is  the  vertical  grid  index.  In  contrast,  the  horizontal 
derivatives  are  solved  using  a  lower-order  upwind  differencing  scheme  (23),  such  that  the  first 
derivatives  are  computed  as 


Wdj) 

dx  Ax 


(34) 


and  the  second  derivatives  are  given  as 

d20(U)  </>(Uj)-2<Ki-\,j)+<ki-2,j) 

dx 2  (Ax)2 

A  set  of  simultaneous  equations  for  each  of  12  computed  variables  is  then  solved  iteratively 
using  the  Thomas  algorithm,  a  tridiagonal  matrix  solver  (24).  In  addition,  the  solution 
implements  a  relaxation  scheme  for  all  the  computed  variables,  similar  to  that  described  by 
Wilson  (22),  i.e., 


(j)  =  a<f>n  +  (l  -  a)</>"  1  ,  (36) 

where  a  =  0.05,  n  is  the  current  iteration,  and  n  -  I  is  the  previous  iteration.  The  relaxation 
scheme  primarily  affects  the  rate  at  which  the  solution  converges,  not  how  well  the  solution 
converges. 

Then,  based  on  the  earlier  works  of  Meyers  and  Paw  U  (12)  and  Katul  and  Albertson  (11),  the 
following  (model)  top  and  bottom  boundary  conditions  are  applied: 

At  z  =  o  : 


(iij  =  (w^  =  0  ;  ^6^  =  303.15  ; 

p)  =  o  ; 

(u =  o  j  (w'd'^=  (ii'0''j  =  o  ; 
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«’2  =  V2  =  W2  =  **  =  o; 


M*  =  (MW 


(friction  velocity  at  z  =  2Az); 


(potential  temperature  scaling  constant  at  z  =  2Az); 


and 

d0) 

ClilU 

k2Az 

dz 

dip) 

0  and  — — - 

=0 . 

dz 

(k  =  0.4  is  von  Karman’s  constant); 


At  z  =  30  m : 


dz  k{z  -d) 


d\d )  0* 
dz  k(z  -d) 

rrVif  (  /7  ~  2  I,  \ 


( d  is  displacement  height  ( d  ~  |  h  )  and  h  is  canopy  height); 


fc  Utop  - 


\  /  (/,/r+lOAz) 
z-d  \ 
h  +  10Az  -  d  ) 


(Utop  =7.0  ms-1  ); 


Qh  =  J  .s>,  dz  (kinematic  heat  flux,  Qh  =  iw'0') ); 


6>»  =  (potential  temperature  scaling  constant); 


5(w)  dfp) 

Ai  =  0  and  AA  =  o; 

&  9z 


u'wj  =  — m*  ;  (yrfO'j  =  Qh  \  and  (u'O'j  =  -3.0  Qh  \ 

u'2\  =  3.5  ;  /v'2\  =  1.5  id  ;  /w'2\  =  1.5  id  ;  and  Iff'2)  =  4.0  0}  . 


In  addition,  zero-gradient  conditions  are  assumed  at  the  lateral  boundaries.  However,  a  weighted 
smoother,  similar  to  that  proposed  by  Mahrer  and  Pielke  (25),  is  applied  to  the  first  three  and  last 
three  horizontal  grid  points  to  reduce  undesired  boundary  effects,  i.e., 

~  0-25(^+ij)  +  l,y')j  •  (32) 

Finally,  lower-order  Newton-Cotes  fonnulas  for  numerical  integration,  i.e.,  the  trapezoidal  rule 
and  Simpson’s  one-third  rule,  are  applied  to  derive  the  wind  flow  streamlines. 
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2.5  Forest  Canopy  Architecture 

As  shown  in  sections  2.2  and  2.3,  canopy  architecture  plays  an  important  role  in  defining  the 
momentum  and  heat  flux  divergence  through  the  forest  layer.  Following  the  discussions  by 
Massman  (26)  and  Meyers  et  al.  (27),  it  was  suggested  that  forest  canopies  may  confonn  to  one 
of  three  general  leaf  area  distribution  profiles,  as  shown  in  Figure  1.  It  is  clear  that  leaf  area 
distributions  are  not  always  symmetric  about  the  layer  of  maximum  foliage  density  (like  profile- 
1)  but  may  be  more  often  skewed  upward  toward  the  top  of  the  forest  canopy.  By  definition,  leaf 

h 

area  index  is  LAI  =  J  A(z)  dz  ,  where  A(z)  is  the  leaf  area  density  through  the  small  vertical  layer 

o 

between  z  and  z  +  Az  per  unit  surface  area  of  ground  below  (28).  Values  for  leaf  area  index  for 
forests  vary  but  have  been  reported  most  often  in  the  range  LAI  =  1  to  5  (29).  In  section  3,  2-D 
model  results  will  be  shown  for  the  case  corresponding  to  profile-2.  Also,  in  this  study  LAI  = 

3.0. 


Figure  1.  Normalized  vertical  profiles  of  leaf  area  density  for  forest 
canopies. 


3.  Model  Results 


3.1  Uniform  Forest  Stands 

In  this  section,  several  initial  model  results  are  presented  for  a  2-D  uniform  forest  stand.  Figure 
2  shows  the  calculated  fields  for  horizontal  wind  velocity,  (uj  ,  in  units  ms  ’,  and  air 

temperature,  (6Sj ,  in  units  °C.  For  this  example,  upper  level  (i.e.,  model  top)  wind  velocity  is 
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umax  =  7.0  ms  ,  leaf  area  index  is  LAI  =  3,  and  forest  canopy  height  is  h  =  10  m.  In  contrast, 
Figure  3  shows  individual  profiles  for  wind  velocity  and  air  temperature  derived  from  the  current 
2-D  calculation  in  comparison  to  profiles  derived  from  an  earlier  1-D  modeling  study  (1,  2). 
Profile  results  from  the  1-D  and  2-D  models  agree  reasonably  well.  Small  differences  are  due 
(possibly)  to  the  expanded  numerical  grid.  Finally,  Figure  4  shows  some  spurious  oscillations 
that  result  when  Av  =  20  m.  Here,  it  is  likely  that  the  solution  is  unstable  because  Ax/Z(N)  <  1, 
where  Z(N)  is  the  height  at  the  model  top.  Based  on  an  analysis  of  the  equation  set  and  several 
additional  calculations  (not  shown),  this  condition  for  stability  appears  to  be  valid.  The 
computed  variable  in  Figure  4  is  horizontal  wind  velocity,  (u)  ,  in  units  ms  '. 
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- 2-D  Model  :  /  I  :l  - 2-D  Model 


0  2  4  6  28  29  30  31 


<  u  >  <  T  > 

_ (a) _ (b) _ 

Figure  3.  Profiles  of  horizontal  wind  velocity  ( (u\  ,  in  units  ms-1)  and 


air  temperature  ( (8j ,  in  units  °C)  derived  from  the  current 

2-D  calculation  shown  in  comparison  to  1-D  profiles  derived 
from  an  earlier  study. 


Figure  4.  Spurious  oscillations  due  to  numerical  instability  in  the 
2-D  model  results  for  uniform  forest  stands.  The 

computed  variable  is  horizontal  wind  velocity,  (ii) ,  in  units 

ms-1. 
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3.2  Nonuniform  Forest  Stands 

Figures  5-7  show  the  current  2-D  model  results  for  a  nonunifonn  forest  stand,  i.e.,  which 
contains  a  single  step  change  in  canopy  height.  The  computed  variables  in  Figure  5  are 

horizontal  wind  velocity,  (uj  ,  in  units  ms-1,  and  the  wind  flow  streamlines.  A  single  step 
change  in  canopy  height  (h)  at  one  grid  increment  past  the  midpoint  (i.e.,  X/2  +  Ax )  is  shown  via 
open  rectangles,  where  h  =  8  m  on  the  left  side  and  h  =  10  m  on  the  right  side.  The  computed 
variables  in  Figure  6  are  vertical  wind  velocity,  (w\,  in  units  ms-1,  and  kinematic  (fluctuation) 

pressure  (p'j ,  in  units  m  V2.  The  computed  variables  in  Figure  7  are  air  temperature,  (O  j ,  in 
units  °C,  and  the  effective  speed  of  sound,  (Cejf)’  in  units  ms-1.  Finally,  Figure  8  shows 

spurious  oscillations  that  result  when  Ax  =  20  m  instead  of  Ax  =  50  m.  Here  again,  a 
computational  instability  is  brought  about  (possibly)  because  A x/Z(N)  <  1 . 


(a)  horizontal  wind  velocity,  (u\  ,  in  units  ms  and  (b)  wind  flow  streamlines.  A  single  step  change  in 

canopy  height  at  X/2  +  Ay  is  shown  using  open  rectangles,  where  h  =  8  m  on  the  left  side  and  h  =  10  m 
on  the  right  side. 
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z(m) 


Figure  6.  Model  results  for  nonuniform  forest  stands:  (a)  vertical  wind  velocity,  no) ,  in  units  ms  1  and  (b) 


kinematic  (fluctuation)  pressure  (j)j  ,  in  units  m  2s  2.  A  single  step  change  in  canopy  height  occurs  at 
X/2  +  Ax,  where  h  =  8  m  on  the  left  side  and  h=  10  m  on  the  right  side. 


Figure  7.  Model  results  for  nonuniform  forest  stands:  (a)  air  temperature,  (8j ,  in  units  °C  and  (b)  effective  speed 

of  sound,  (Cejj  ^ ,  in  units  ms-1.  A  single  step  change  in  canopy  height  occurs  at  X/2  +  Ax,  where  h  =  8  m 
on  the  left  side  and  h  =  10  m  on  the  right  side. 
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Figure  8.  Spurious  oscillations  due  to  numerical  instability  in  the 
2-D  model  results  for  nonuniform  forest  stands.  The 

computed  variable  is  horizontal  wind  velocity,  (uj  ,  in 

units  ms-1. 


Figures  9-11  show  the  current  2-D  model  results  for  a  second  nonunifonn  forest  stand.  In  this 
example,  a  larger  step  in  canopy  height  is  incorporated  at  the  lower  boundary,  where  h  =  4  m  on 
the  left  side  and  h=  10  m  on  the  right  side.  The  computed  variables  in  Figure  9  are  horizontal 

wind  velocity,  (uj  ,  in  units  ms”1,  and  the  wind  flow  streamlines.  The  computed  variables  in 
Figure  10  are  vertical  wind  velocity,  (w^j ,  in  units  ms”1,  and  kinematic  (fluctuation)  pressure  (jjj, 
in  units  m  2s”2.  The  computed  variables  in  Figure  1 1  are  air  temperature,  (Oj ,  in  units  °C,  and 


the  effective  speed  of  sound,  (Cejj \,  in  units  ms  '.  Finally,  Figure  12  shows  that  care  needs  to  be 


taken  when  applying  numerical  (integration)  schemes  across  sharp  discontinuities  (e.g.,  to 
calculate  the  streamlines).  Here,  the  use  of  Simpson’s  one-third  rule  for  both  integrals,  i.e., 


j d i//  =  - j  wdx  +  J udz  ,  resulted  in  a  computational  instability  (i.e.,  2  Ax  numerical  waves).  Upon 


further  analysis,  these  were  found  to  be  contained  mainly  in  the  solution  to  -  wdx .  A  stable 


solution  was  later  obtained  by  applying  a  simpler  trapezoidal  rule  to  solve  this  integral.  Note 
that  while  Ax  =  50  m  in  Figure  12,  the  number  of  horizontal  grid  points  is  expanded  (i.e.,  MPT 


=  40). 
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z(m) 


Figure  9.  Model  results  for  nonuniform  forest  stands,  i.e.,  those  that  contain  a  single  step  change  in  canopy  height: 

(a)  horizontal  wind  velocity,  ,  in  units  ms-1,  and  (b)  wind  flow  streamlines.  A  single  step  change  in 

canopy  height  at  X/2  +  Ax  is  shown  using  open  rectangles,  where  h  =  4  m  on  the  left  side  and  h  =  10  m  on 
the  right  side. 


Figure  10.  Model  results  for  nonuniform  forest  stands:  (a)  vertical  wind  velocity,  liv\ ,  in  units  ms  and  (b) 


kinematic  (fluctuation)  pressure  (jjj  ,  in  units  m  2s  2.  A  single  step  change  in  canopy  height  occurs  at 
X/2  +  Ax,  where  h  =  4  m  on  the  left  side  and  h  =  10  m  on  the  right  side. 
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;ure  11.  Model  results  for  nonuniform  forest  stands:  (a)  air  temperature,  (@j ,  in  units  °C,  and  (b)  effective 

speed  of  sound,  {Ceg\,  in  units  ms-1.  A  single  step  change  in  canopy  height  occurs  at  X/2  +  Ax,  where 
/?  =  4mon  the  left  side  and  h  =  10  m  on  the  right  side. 


Figure  12.  Numerical  instability,  i.e.,  2Ax  waves,  in  the  computation  of  wind  flow  streamlines  for  nonuniform 

forest  stands:  (a)  J  d  (/  =  -J  wdx  +  J  udz  and  (b)  Jrf  y/  =  -j" wdx  only.  A  single  step  change  in  canopy 
height  occurs  at  X/2  +  Ax,  where  h  =  4  m  on  the  left  side  and  h  =  10  m  on  the  right  side. 
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4.  Summary  and  Conclusions 


A  new  and  relatively  unique,  physics-based,  meteorological  computer  model  for  2-D  forest 
canopy  wind  flow,  temperature,  and  turbulence  calculations  has  been  developed.  The  current 
2-D  model  is  based  on  the  same  basic  conservation-law  equations  and  (second-order  turbulence 
closure)  modeling  assumptions  that  were  implemented  in  an  earlier  1-D  model  study  to 
mathematically  represent  the  mechanical  and  thennodynamic  influences  on  the  speed  of  sound  in 
the  forest  environment.  The  2-D  computer  model  has  been  implemented  in  Fortran.  Valid 
numerical  techniques  from  earlier  works  have  been  applied  to  solve  for  each  of  twelve  computed 
variables,  to  include  the  mean  flow  vertical  velocity  and  the  kinematic  (fluctuation)  pressure. 

The  horizontal  derivatives  were  solved  using  lower-order,  upwind  differencing.  The  vertical 
derivatives  were  solved  using  lower-order,  central  differencing.  Second-order,  ordinary 
differential  equations  for  the  profiles  were  solved  using  a  tridiagonal  matrix  algorithm.  Thus, 
several  satisfactory  solutions  were  achieved  for  uniform  and  nonuniform  forest  stands  (for  coarse 
horizontal  grid  spacing,  i.e.,  Ax  =  50  m).  In  addition,  a  valid  condition  for  numerical  stability 
was  determined,  as  A x/Z(N)  >  1 . 

In  future  modeling  works,  we  will  continue  to  investigate  these  and  alternate  numerical  schemes, 
which  are  accurate,  fast,  and  robust,  to  solve  the  computed  fields  within  and  above  a  realistic 
forest  canopies,  particularly  those  containing  sharp  discontinuities  at  the  lower  boundary  (e.g., 
forest  edges  [JO]).  Also,  the  presence  of  hills  may  significantly  alter  the  flow  field  inside 
canopies  (31).  Most  turbulence  models,  including  the  model  described  herein,  are  for  flat 
surfaces.  Therefore,  additional  2-D  forest  canopy  models  may  begin  to  consider  uneven  terrain. 
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